Critical Properties of the Kitaev-Heisenberg Model 
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We study the critical properties of the Kitaev-Heisenberg (KH) model on the honeycomb lattice at 
finite temperatures that might describe the physics of the quasi two-dimensional (2D) compounds, 
Na2lr03 and Li2lr03. The model undergoes two phase transitions as a function of temperature. At 
low temperature, thermal fluctuations induce magnetic long-range order by the order-by-disorder 
mechanism. This magnetically ordered state with a spontaneously broken Z§ symmetry persists 
up to a certain critical temperature. We find that there is an intermediate phase between the 
low-temperature, ordered phase and the high-temperature, disordered phase. Finite-sized scaling 
analysis suggests that the intermediate phase is a critical Kosterlitz-Thouless (KT) phase with 
continuously variable exponents. We argue that the intermediate phase has been observed above 
the low-temperature, magnetically ordered phase in Na2lrC>3, and also likely exists in Li2lr03. 



Introduction. The Ir-based transition metal oxides, in 
which the orbital degeneracy is accompanied by a strong 
relativistic spin-orbit coupling (SOC), have recently at- 
tracted a lot of theoretical and experimental attention 
[ll-Q. This is because the strong SOC creates a different, 
and frequently novel, set of magnetic and orbital states 
due to the unusual anisotropic exchange interactions be- 
tween localized moments which are in turn determined 
by the combination of spin and lattice symmetries. The 
spin-orbital models that describe the low-energy physics 
of iridium systems often include anisotropic terms that do 
not reduce to the conventional easy-plane and easy-axis 
anisotropies because they involve the products of differ- 
ent components of multiple spin operators. These terms 
are responsible for exotic Mott-insulating states , topo- 
logical insulators U | , spin-orbital liquid states 0,13 > 
and non-trivial long-range magnetic orders 0, 0, EJ . 

A prominent example of such an anisotropic spin- 
orbital model is the KH model on the honeycomb lattice 
12l 13 1 which likely describes the low-energy physics of 
the quasi 2D compounds, Na 2 Ir03 and Li 2 Ir03. In these 
compounds, Ir 4+ ions are in a low spin 5d 5 configura- 
tion and form weakly coupled hexagonal layers [1, |g, Q . 
Due to strong SOC, the atomic ground state is a dou- 
blet where the spin and orbital angular momenta of Ir + 
ions are coupled into J c g — 1/2. It was suggested [12l.ll3| 
that the interactions between these effective moments can 
be described by a spin Hamiltonian containing two com- 
peting nearest neighbor (NN) interactions: an isotropic 
antiferromagnetic (AF) Heisenberg exchange interaction 
and a highly anisotropic ferromagnetic (FM) Kitaev ex- 
change interaction [14| . This competition can be de- 
scribed with the parameter, < a < 1, which sets the 
relative strength of these two interactions. At a = 0, the 
coupling corresponds to the AF Heisenberg interaction, 
and at a — 1, it corresponds to the Kitaev interaction. 

This model immediately attracted a lot of attention; 
several theoretical studies were published in the last few 
years [l3|, [T^- 12 1 on both the ground state and its prop- 
erties at finite temperature. The ground state phase di- 



agram of the KH model exhibits three distinct phases: 
the AF Neel phase for small a £ (0., 0.4), the stripy AF 
phase for intermediate a £ (0.4, 0.86), and the disordered 
spin-liquid phase at large a £ (0.86, 1.). While the phase 
transition between the Neel and the stripy phase appears 
to be discontinuous, numerical studies including density 
matrix renormalization group |15| and exact diagonaliza- 
tion results [l3| suggest that the transition between the 
spin liquid and the stripy state is continuous or weakly 
first-order. Additionally, quantum fluctuations select all 
of the magnetically ordered phases to have the order pa- 
rameter point along one of the cubic axes. 

In this Letter, we discuss finite temperature properties 
of the KH model on the honeycomb lattice. A first step 
in this direction was made in Ref . [lj| , where the critical 
ordering scale for the magnetically ordered states was an- 
alyzed using a pseudofermion functional renormalization 
group approach. Here we present numerical results ob- 
tained using Monte Carlo (MC) simulations. We study 
the classical KH model because the corresponding quan- 
tum model has a sign problem precluding quantum MC 
analysis and also because the existence of long-range or- 
der at low temperatures in Na2±r03 and in Li2ir03 in- 
dicates that quantum fluctuations are not dominating in 
these materials [HI]. 

We show that the thermal fluctuations of classical spins 
give rise to two distinct temperature dependent effects. 
At low temperature they predominantly act as the source 
of the order-by-disorder phenomenon and select collinear 
magnetic order where the spins are oriented along one 
of the cubic directions. There are six possible ordered 
states, one of which is spontaneously chosen by the sys- 
tem. At high temperatures, when T is larger than any 
energy scale in the system, the fluctuations destroy any 
order putting the KH model into a three dimensional 
paramagnetic state. The main goal of our study is to see 
how these two phases are connected. 

We argue that the classical KH model effectively be- 
haves like a six-state clock model 18-211 and that it un- 
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FIG. 1: Four possible magnetic configurations: (a) the FM 
ordering; (b) the two-sublattice, AF Neel order; (c) the stripy 
order; (d) the zigzag order. Open and filled circles correspond 
to up and down directions of spins. 
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FIG. 2: Histograms of the order parameter mjv(s)i obtained 
for the system with 2*84*84 spins in the ordered phase, (a) 
and (e), in the intermediate phase, (b)-(c) and (f)-(g), and 
in the disordered phase, (d) and (h). Histograms (a)-(d) are 
computed for a — 0.25, and (e)-(h) are for a = 0.75. The 
histograms are presented on the complex plane (Re |m]v(s)|, 
Im \m N(s) \). 



of temperature separating three phases: a low-T ordered 
phase, an intermediate critical phase, and a high-T dis- 
ordered phase. The critical phase has an emergent, con- 
tinuous U(l) symmetry which is fully analogous to the 
low-T phase of the XY model, a well-known KT phase 
of critical points with floating exponents and algebraic 
correlations. Here we present numerical data only for 
a = 0.25 and a = 0.75 since these values likely character- 
ize the ratio between the AF Heisenberg interaction and 
the Kitaev interaction in Na2lr03 and Li2lr03. However, 
we note that recent inelastic neutron scattering measure- 
ments on Na2lr03 have shown that the KH model alone 
is insufficient to describe the magnetic properties of this 
compound It lias been demonstrated that it is es- 
sential to include substantial further-neighbor exchanges 
to describe both the zigzag ground state and the exci- 
tation spectrum in Na2lr03. The full finite-temperature 
phase diagram for the KH model with second and third 
neighbor exchange interactions will be published else- 
where m. 

The Model. The classical version of the KH model 
which describes the interactions among the J = 1/2 de- 
grees of freedom of Ir + ions reads as 

H = -J K ^2S7SJ + J H ^S i S J . (1) 

where the spin quantization axes are taken along the cu- 
bic axes of the Ir06 octahedra. 7 = x,y, z denotes the 
three bonds of the honeycomb lattice. The exchange con- 
stants, Jk = 2a and Jh = 1 — a, correspond to the Ki- 
taev and Heisenberg interactions which can be derived 
from a multiorbital Hubbard Hamiltonian [l3j |. 

Order by Disorder. The symmetry of the KH model 
combines the cubic symmetry of both the spin and the 
lattice space. It consists of simultaneous permutations 
between the x, y, z spin components and a C^-rotation of 
the lattice which defines a discrete symmetry. The classi- 



cal ground state has a higher symmetry than that of the 
Hamiltonian - the ground state energy does not change 
under a simultaneous rotation of all spins. Since this 
applies only to the ground state, the KH model has only 
an accidental continuous rotation symmetry. Its actual 
symmetry is discrete; at zero temperature, the "pseudo" 
SU(2) symmetry is broken by quantum fluctuations that 
restore the underlying cubic symmetry of the model 
The magnetically ordered phase is gapped with a spin gap 
that corresponds to the finite energy cost of deviating the 
order parameter from one of the cubic axes. We show in 
the following that thermal fluctuations of classical spins 
at finite T also select a collinear spin configuration whose 
order parameter points along one of the cubic axes. 

Parameters of the Simulations. We have carried out 
classical MC simulations of the model (TIJ using the 
standard Metropolis algorithm. In our MC simulations, 
we treat the spins as three-dimensional (3D) vectors, 
Si = (Sf , Sf , Sf), of unit magnitude with (Sf) 2 + (S?) 2 + 
(Sf) 2 = 1 at every site. At each temperature, more than 
10 7 MC sweeps were performed. Of these, 5 * 10 5 were 
used to equilibrate the system, and afterwards only 1 out 
of every 5 sweeps was used to calculate the averages of 
physical quantities. We present all energies in the units 
of Jh and assume ks = 1. The calculations were carried 
out on several finite systems with size 2 * L * L that are 
spanned by the primitive vectors of a triangular lattice 
ai = (1/2, v3/2) and &2 — (1,0) with a 2-point basis 
using periodic boundary conditions. 

Results. To study the possible phases of the model HJ, 
we introduce four magnetic configurations (Fig.[T]): a FM 
order, a simple two-sublattice AF Neel order, a stripy or- 
der, and a zigzag spin order. The classical energies of 
these states can be easily computed: E^ 1 = 3 — 5a, 
E$ = -3a + 1, E% = -a - 1, and E% = 5a - 3 for 
the FM, the zigzag, the stripy and the Neel phases, re- 
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FIG. 3: The log- log plots of the order parameter mjv(s) 
as a function of system size L at various temperatures. The 
solid curves indicate the linear behavior that corresponds to a 
power-law dependence, mjv(s) ~ L~ v ^ 2 , cooresponding to the 
intermediate critical phase. The dashed curves show deviation 
away from the linear behavior outside the critical phase. 



spectively. For < a < 1, the classical ground state 
is either the Neel AF with the vector order parameter 
Af = — S^) or the stripy phase described 

by $ = E i= „(Su - S lB + S lC - S iD ). Here, A, B 
and A, B, C, D denote either two or four sublattices that 
respectively characterize the Neel AF and stripy order. 
The classical phase transition between them occurs at 
a = 1/3. At a = 1, the FM, stripy, and zigzag phases 
all have the same classical energy. However, the classi- 
cal degeneracy of this point, which corresponds to the 
pure Kitaev model, is much higher. This limit has been 
thoroughly studied by Baskaran et al. (22^ . 

To make an analogy to the six-state clock model, we 
map the order parameter describing the magnetically or- 
dered phase of the KH model onto a 2D complex order 
parameter, m^s) = Yli=i \ m i,N(s) such that the six 
possible ordered states are characterized by 0$ = irrii/3, 
rii = 0,..5 20]. The mapping is exact only well within 
the ordered state since there is no guarantee that the 
thermal fluctuations of the order parameter will actually 
have a 2D character given that the spin degrees of free- 
dom are three-dimensional. Depending on the strength 
of the spin stiffness in different directions, the long-range 
low-T magnetic order can be destroyed in one of several 
ways. If the stiffness of thermal fluctuations along the 
circle is softer than the stiffness of fluctuations in the di- 
rection transverse to the circle, the long-range order may 
be destroyed by a discontinuous first-order transition, by 
two continuous phase transitions with an intermediate 




FIG. 4: A snapshot of the coarse-grained order parameter 
(mjv) at T = 0.168. The vortex-like topological excitations 
are evident. 



partially ordered phase, or by two K T p hase transitions 
with an intermediate critical phase [18l42l| . In the last 
scenario, the critical phase is destroyed by topological ex- 
citations in the form of discrete vortices whose existence 
is directly related to the emergence of a continuous sym- 
metry; the high-T transition will first bring the system 
into a disordered phase where fluctuations are primarily 
2D, and the crossover to the 3D paramagnet occurs at 
even higher temperatures. 

In Fig. [2] we present the results of the histogram 
method for the complex order parameter. At low temper- 
atures, Figs. [2] (a) and (e), a sixfold degeneracy present 
in the ordered phase is seen. For both a = 0.25 and 
a = 0.75, the six states which have the highest weight 
in the histogram are where the order parameter mjv(s) 
points along one of the cubic axes. In Figs. [2] (b) and (f), 
when the temperature increases beyond a certain critical 
temperature, a continuous t/(l) symmetry emerges sig- 
naling both the disappearance of the sixfold anisotropy 
and the appearance of the critical phase. The forma- 
tion of vortices can be seen in Fig. [4] where we present a 
snapshot of the coarse-grained order parameter (ro/y) at 
T = 0.168. Upon a further increase in temperature, the 
amplitude of the order parameter decreases (Figs. [2] (c) 
and (g)) until it shrinks to zero indicating the transition 
to the paramagnetic phase (Figs. [2] (d) and (h)). 

To better understand the properties of the intermedi- 
ate phase and to confirm its critical nature, we performed 
the finite-size scaling analysis appropriate for KT tran- 
sitions 



241 ] . The full finite-size scaling analysis is rather 
involved and will be reported elsewhere [23J. Here we 
present only the scaling behavior of the order parame- 
ter. At the KT transition, the order parameter exhibits 
the power law dependence on system size, m ~ L~''/ 2 . 
As each point of the intermediate critical phase can be 
understood as a critical point, the power law behavior 
of the order parameter should hold throughout the en- 
tire phase. We found that the boundaries of the critical 
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FIG. 5: The Binder cumulant as a function of temperature 
for (a) a = 0.25 and (b) a — 0.75. From the crossing points of 
different Binder's curves, we estimate T C1 — 0.152 and T C1 = 
0.124 for a = 0.25 and a — 0.75, respectively. 



phase are characterized by critical exponents close to 1 /9 
and 1/4 for the lower and upper boundaries at T Cl and 
T C2 , which is in agreement with critical exponents for 
the six-state clock model obtained by the renormaliza- 
tion group analysis [lij]. Fig. [3] shows the log- log plots of 
the order parameter m^s) as a function of system size 
for different temperatures. For a — 0.25, the data points 
in Fig. [3] a) show a linear behavior in the temperature 
interval between T Cl ~ 0.152 and T C2 ~ 0.162, in which 
there are several critical lines characterized by r) between 
1/9 and 1/4. For a = 0.75, we have detected the critical 
phase in the temperature interval between T C1 ~ 0.125 
and T C2 ~ 0.127. 

The lower transition temperature T Cl can be indepen- 
dently determined using fourth-order Binder cumulant 
(Figs. 0(a) and (b)). The Binder cumulant has a scaling 
dimension of zero; thus the crossing point of the cumu- 
lants for different lattice sizes provides a reliable estimate 
for the value of the critical temperature T Cl at which the 
long range order is destroyed. The crossing points for 
a = 0.25 and a = 0.75 are T Cl = 0.152 and T Cl = 0.124, 
respectively. They are in good agreement with estimates 
obtained from the log-log plots in Fig. [3] 

In Figs. M (a) and (b) we present the temperature de- 
pendence of the specific heat, C = ((E 2 ) - (E) 2 )/NT 2 . 
While the low-T transition, seen as small peak at tem- 
peratures T C1 = 0.152 and 0.1247 for a = 0.25 and 0.75, 
respectively, is in a good agreement with our previous 
estimates, the features corresponding to the high-T tran- 
sition T C2 are barely distinguished by eye. This is not 
surprising as the high-T transition is a usual KT tran- 
sition at which the specific heat does not diverge at the 
critical point |26|. It is also likely that the high-T KT 
transition might be shadowed by the crossover to the 3D 
paramagnet, which is seen in Fig. [5] as a very broad hump 
at higher-T. 

Our findings for the specific heat show a lot of simi- 
larities between the experimental data obtained on the 
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FIG. 6: Specific heat C as a function of temperature for (a) 
a = 0.25 and (b) a = 0.75. 



Na2lrC>3 and Li 2 Ir03 compounds by Refs. @, HI an d @|- 
In Na2lrC>3, both the lambda- like anomaly at the Neel or- 
dering temperature, Tn = 15 K, and a broad tail which 
extends into higher temperatures are seen in the specific 
heat measurements [4j. The latter suggests the presence 
of short-range order above the bulk 3D ordering that can 
be understood by our proposed scenario of the critical 
phase. 

Let us estimate the temperatures of the KT transitions 
and the width of the critical phase in Na2lrC>3 based on 
our results obtained for the KH model with a = 0.25. On 
the mean field level, the exchange on the NN bonds may 
be estimated from the classical energy, J\ ~ (3 — 5a)/3, 
in the Neel phase. From the recent neutron scattering 
experiment [7|], the NN exchange in Na2lr03 was esti- 
mated to be Ji = 4.17 meV. In the bulk of our paper, 
all energies were measured in the units of Jh, and thus 
we estimate J\ to be equal to 12.7 meV. This gives the 
prediction for the critical temperature to be T C1 = 16.8 
K, which is very close to the experimental value Tjy = 15 
K [3, [1]. Our estimate for the upper boundary of the 
critical phase is T C2 = 17.7 K which makes the predicted 
critical phase very narrow. We note here that the critical 
phase survives in the extended KH model with included 
further-neighbor exchange couplings [B|, 0, l25| which are 
essential for comparison with experiment. However, in 
order to determine the upper boundary of the critical 
phase additional extensive numerical simulations must 
be performed. 
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